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PACS 47.52. + j - Chaos in fluid dynamics 

PACS 47.20. Bp - Buoyancy-driven instabilities (e.g., Rayleigh-Benard) 

PACS 47.27. ed - Dynamical systems approaches 

Abstract -Transition to hyperchaotic regimes in Rayleigh-Benard convection in a square period¬ 
icity cell is studied by three-dimensional numerical simulations. By fixing the Prandtl number at 
P = 0.3 and varying the Rayleigh number as a control parameter, a bifurcation diagram is con¬ 
structed where a route to hyperchaos involving quasiperiodic regimes with two and three incom¬ 
mensurate frequencies, multistability, chaotic intermittent attractors and a sequence of boundary 
and interior crises is shown. The three largest Lyapunov exponents exhibit a linear scaling with 
the Rayleigh number and are positive in the final hyperchaotic attractor. Thus, a route to weak 
turbulence is found. 
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c/D Introduction. — Thermal convection refers to the 
^motion of a heat-conducting fluid due to the presence of 


oi 


temperature differences. Convective flows are of interest 
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.n many areas ranging from technological processes (cool¬ 
ing of electronic devices, drying, material processing) to 
natural phenomena (convection in the terrestrial atmo- 
sphere, oceans, mantle, outer core and stellar convection). 
(^Origin of magnetic fields of many astrophysical objects 
80 is explained in the framework of the dynamo theory [l], 
O* where magnetic fields are generated by convective motions 
*^of electrically conducting fluids in their interior; chaotic 
80 convective flows play an important role in the fast dynamo 
< ^ > theory. 

Thermal convection in the dimensionless form is char- 
• • acterised by the Rayleigh number (Ra), that measures the 
magnitude of thermal buoyancy force, and the Prandtl 
number (P), the ratio of kinematic viscosity to thermal 
Vh diffusivity. In the context of Astrophysics, values of the 
^Prandtl number vary from small (for the solar convec¬ 
tive zone, P ~ 10“ 7 [ 2 ]) to large (for mantle convection, 
P ~ 10 23 [ 3 ]) values. The Prandtl number for the terres¬ 
trial outer core has intermediate values estimated to be 
between 0.1 and 0.5 [4][5], however, in many terrestrial 
convective dynamo simulations a larger value (P ~ 1) is 
employed (see, e.g., [6] and references therein) using the 
argument that it expresses the “effective” or ” turbulent” 
value of the diffusivity [7]. 

In the dynamo theory, interest in the use of realistic val¬ 


ues for the Prandtl number was enhanced by the following 
findings. Varying the value of the Prandtl number from 0.1 
to 10 was found to have a strong influence on the morphol¬ 
ogy and dynamics of convection in the Earth’s outer core 
[8]. Strong dependence of the magnetic fields generated 
by convective flows on the value of the Prandtl number 
was found for 0.2 < P < 5 [9 . In [10 , moderately low 
Prandtl numbers were beneficial for magnetic field gen¬ 
eration in rotating spherical shells and special attention 
was devoted to dynamos for P = 0.1. Analysis of the 
kinematic dynamo problem in [ll] showed that convective 
attractors for P = 0.3 are beneficial for the magnetic field 
generation (in comparison to the attractors for P = 1 and 
6.8). 

Most of the convective flows in nature are turbulent. 
Thermal convection in a plane horizontal layer, called 
Rayleigh-Benard convection, has been used for decades 
as one of the simplest examples of realistic hydrodynamic 
systems driven out of equilibrium where the simple flow 
becomes complex (turbulent) in a variety of bifurcation se¬ 
quences, revealing different mechanisms of instability and 
demonstrating many common nonlinear dynamics phe¬ 
nomena, e.g., spontaneous symmetry breaking, pattern 
formation, intermittency, synchronisation, etc. (for a re¬ 
view, see [T2l13 ). A common approach to study transition 


to turbulence in the framework of the dynamical systems 
theory is to study the evolution from simple non-chaotic 
(steady, time-periodic and quasiperiodic) attractors of the 
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convective system to the chaotic ones. In chaotic attrac¬ 
tors, trajectories are sensitive to initial conditions, i.e., 
initially close trajectories diverge in time, which is quan¬ 
titatively characterised by the Lyapunov spectrum. Tur¬ 
bulent systems usually display hyperchaos, i.e., more than 
one positive Lyapunov exponent. 

There is a wide range of works on bifurcation analysis 
of Rayleigh-Benard convection as a function of Ra and P, 
most of which are devoted to the formation and destabil¬ 
isation of steady planar convective rolls [Id 17 and tran¬ 
sition to chaos 
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In 


a bifurcation diagram 


is presented where a series of steady states representing 
different patterns is numerically obtained as a function 
of Ra in a small cylindrical domain. Also in cylindri¬ 
cal domains, hyperchaotic states were found in studies of 
spiral defect chaos using simulations of three-dimensional 
Rayleigh-Benard convection, where the spectrum of Lya¬ 
punov exponents was used to quantify extensivity in spa- 


tiotemporal chaos 20 


Still regarding transition to hyperchaos in Rayleigh- 
Benard convection, in 22 , a square periodicity cell with 
aspect ratio L = 2y/2 and stress-free boundaries was stud¬ 
ied for P = 6.8 in two-dimensional convection for Ra 
up to 32875 using a low-dimensional model (16 Fourier 
modes); the sequence of regimes is as follows: time- 
periodic, quasiperiodic with two basic frequencies, phase- 
locked (periodic) and then chaotic state (some of the 
chaotic states are hyperchaotic with two positive Lya¬ 
punov exponents). 

In this letter we study the transition to hyperchaotic 
regimes (also referred to as transition to weak turbulence 
[23] ) in three-dimensional Rayleigh-Benard convection for 
P = 0.3 in a square convective cell for Ra increasing from 
1720 (time-periodic state) to 2500 (hyperchaotic state). 
We show that a sequence of crises involving quasiperiodic 
and chaotic attractors, as well as chaotic saddles, is re¬ 
sponsible for the evolution of the attractors of the system 
from periodic convective states to hyperchaos with at least 
three positive Lyapunov exponents. 


Statement of the problem and solution. — A 

Newtonian fluid flow in a horizontal plane layer is con¬ 
sidered, where the fluid is uniformly heated from below 
and cooled from above. Fluid flow is buoyancy-driven and 
the Boussinesq approximation is assumed. We adopt the 
vertical size of the layer as a length scale, the vertical heat 
diffusion time as a time scale, and the vertical temperature 
gradient as a temperature scale. Then, in a Cartesian ref¬ 
erence frame with the orthonormal basis (ei, e 2 , 63 ), where 
e 3 is opposite to the direction of gravity, the equations 
governing the convective system are (see, e.g., p~5]): 

dv/dt = PV 2 v + v x (V x v) + PRa#e 3 — Vp, (1) 
99/dt = V 2 <9 - (v • V)0 + w 3 , ( 2 ) 

V • v = 0, (3) 

where v(x,t) = (tq, v 2 , v 3 ) is the fluid velocity, p(x, t) the 

pressure, and 0(x, t) = T(x, t) — (Tj. + (T 2 — Ti)xs) is the 


difference between the temperature T and its linear profile. 
Temperatures of the horizontal boundaries at the bottom, 
Ti, and top, T 2 , are maintained constant, with T\ > T 2 . 
Here, the spatial coordinates are x = (aq, x 2 , £ 3 ) and t 
stands for time. The non-dimensional parameters are the 
Prandtl number, P, and the Rayleigh number, Ra. 

The horizontal boundaries of the layer are assumed to be 
stress-free, dvi/dxs = dv 2 /dx% = V 3 = 0 , and maintained 
at constant temperatures, 6 = 0. A square convective cell 
is considered, x £ [0, L] 2 x [0,1], and all the fields are pe¬ 
riodic in the horizontal directions, aq and a; 2 , with period 
L. The linear theory 114] suggests that at the onset of 
convection in an infinite layer (Ra c =657.5 for the bound¬ 
ary conditions under consideration) the critical horizontal 
wavenumber is tt/V^, independently of the Prandtl num¬ 
ber; here, L = 4 is taken, hence the most unstable mode 
at the onset is aligned with the diagonal of the cell. 

Following [II], we studied attractors of the convective 
system for P = 0.3 and 1720 < Ra < 2500. The attrac¬ 
tors were obtained integrating the system in time starting 
from an attractor for a neighboring value of Ra (in most 
cases located at distance 10); the first 1500 eddy turnover 
times of the largest eddies were disregarded as transients. 
To check if multiple attractors co-exist, for several values 
of Ra the problem was solved for four initial conditions 
defined by random Fourier coefficients of v and 0 with ex¬ 
ponentially decaying spectrum and the following values of 
kinetic energy: E v { 0) = 1, 100, 400, 2500. 

For a given initial condition, the system under consid¬ 
eration is integrated numerically forward in time using 
the standard pseudospectral method [24]: the fields are 
represented as Fourier series in all spatial variables (ex¬ 
ponentials in the horizontal directions, sine/cosine in the 
vertical direction), derivatives are computed in the Fourier 
space, multiplications are performed in the physical space, 
and Orszag’s 2/3-rule is applied for dealiasing. The sys¬ 
tem of ordinary differential equations for the Fourier co¬ 
efficients is solved using the third-order exponential time- 
differencing method ETDRK3 [25] with constant time step 
At = 5 • 10 -4 . The spatial resolution was chosen to be 
32 x 32 x 16 Fourier harmonics (multiplications were per¬ 
formed on a uniform 48 x 48 x 24 grid). For all solutions, 
the time-averaged energy spectra of the velocity decay at 
least by 5 orders of magnitude. For each branch of at¬ 
tractors, several runs with doubled spatial and temporal 
resolution showed no qualitative difference. 

The three largest Lyapunov exponents, Ai > A 2 > A 3 , 
were computed using the technique described in 126], with 
recourse to the operator of linearisation of the governing 
equations 0-0 and the Gramm-Shmidt orthonormali¬ 
sation. Translational invariance of the convective system 
in the horizontal directions gives rise to two vanishing Lya¬ 
punov exponents, which were disregarded in computations 
by removing the corresponding components of the pertur¬ 
bations (note that for non-steady attractors at least one 
vanishing Lyapunov exponent remains since eqs. 0-0 
are an autonomous system). According to the general the- 


p-2 


























Route to hyperchaos in Rayleigh-Benard convection 




0 100 200 1 300 400 50( 0 0.5 1.0 


Fig. 2: Kinetic energy evolution and the Poincare section for 
the quasiperiodic attractor at Ra=2075 from the family A 3 . 
Left panel: kinetic energy E v versus time; right panel: Poincare 
section defined by = 0.25 on the |vi,i,i|) quad¬ 

rant. 


Fig. 1: Time-averaged kinetic energy versus Rayleigh number 
of the attractor families. IC denotes an interior crisis, H a Hopf 
bifurcation. 


ory 27 , for a stable limit cycle (time-periodic attractor) 


Ai = 0 , A 2 , A 3 < 0 ; for a &-torus (quasi-periodic attrac¬ 
tor with k basic frequencies) the largest k Lyapunov ex¬ 
ponents vanish; chaotic and hyperchaotic attractors are 
characterised by at least one and two positive Lyapunov 
exponents, respectively. 

In our computations we traced the kinetic energy, 
E v (i) = f 0 L f 0 L Si v 2 dx/(2L 2 ), as well as the Fourier har¬ 
monics of the flow velocity, Vk it) = f° r some 

wave vectors, k = (k±, & 2 , ^ 3 ). Poincare sections were 
constructed for the Fourier harmonic of the fluid velocity 
for k = ( 1 , 1 , 1 ) on the quadrant (|?5^|, where inter¬ 
section of the trajectory with the plane |uj| = 0.25 was 
considered. Using the solenoidality condition § for the 
Fourier coefficient ^ 1 , 1,1 and the inequalities \zi\ — \z 2 \ < 
\zi + Z 2 I < \zi\ + |^ 21 5 Vzi,Z 2 £ C, one proves that all the 
points on the Poincare section belong to the semi-infinite 
strip |l^i 11 1 — \v\ 11 11 < 0.5, \vl 11 | + |vf 11 | > 0.5. Ab¬ 
solute values eliminate drifting frequencies from consider¬ 
ation. On the Poincare section, a time-periodic attractor 
and a 2 -torus appear as a finite set of points and a curve, 
respectively. 

Results. — In what follows, attractors of the convec¬ 
tive system for 1720 < Ra < 2500 are discussed. This 
range of Ra was also considered in 11 , where bifurca¬ 


tions of the convective attractors were studied for 657.5 < 
Ra < 2500, with transition to chaotic attractors following 
the sequence: periodic-quasiperiodic-chaotic. However, no 
detailed study of the chaotic attractors was performed. In 
the present letter, results of a detailed analysis of the tran¬ 
sition to chaos are reported. Branches of attractors found 
for different intervals of Ra are shown in fig. [lj 

Attractors of the branch Ai are periodic states, exist¬ 
ing for 1720 < Ra < 1780, and are symmetric at any 
time with respect to rotation about the vertical axis by 
7r/2. On increasing Ra, the symmetry is broken and the 
branch A 2 of time-periodic states (with trivial symmetry 
group) drifting along the horizontal directions emanates. 
Although the convective regimes constituting A 2 are for¬ 
mally quasiperiodic, for simplicity, we classify them as pe¬ 


riodic (also referred to as relative periodic orbits [28]) since 
they are periodic in a co-moving reference frame. In what 
follows, we ignore drift frequencies when classifying an at¬ 
tractor as periodic or quasiperiodic. 

A branch of quasiperiodic attractors, A 3 , exists for 
2075 < Ra < 2130, with three basic frequencies (the 
Poincare plane shown on fig. [2] confirms that the regime 
has at least three incommensurate frequencies). We have 
also checked that for all attractors from this branch all the 
three largest Lyapunov exponents vanish. 

Coexisting with A 2 , there is a family of quasiperiodic 
and chaotic attractors denoted by Bi which, on increasing 
Ra, gains stability at Ra=1880. For 1880 < Ra < 1900 the 
attractors in this family are quasiperiodic with two basic 
(incommensurate) frequencies (see fig. [3] (a)). On increas¬ 
ing Ra a period doubling bifurcation occurs, cf. fig. [3] 
(a) and (b), whereby the lowest frequency is halved; this 
quasiperiodic attractor exists in 1910 < Ra < 1980. As Ra 
is increased further, a sequence of chaotic and quasiperi¬ 
odic regimes is observed: for Ra=1990 the convective at¬ 
tractor is chaotic (fig. [ 3 ] (c)); for Ra=2000 it is quasiperi¬ 
odic (fig. [3] (d)); for 2010 < Ra < 2040 regimes are chaotic; 
for Ra=2050 they are quasiperiodic (fig. [3] (h)) and finally, 
for 2060 < Ra < 2070 they are chaotic (fig. [ 3 ] (i) and (j)). 
Measurements of the three largest Lyapunov exponents re¬ 
veal that all chaotic attractors in this sequence have one 
positive and two small in modulus Lyapunov exponents. 
Although the identification of each bifurcation occurring 
in this interval is out of the scope of our study, we note 
that i) windows of quasiperiodicity can be attributed to 
frequency locking [27], ii) some chaotic attractors in the 
sequence display intermittent behaviour, i.e., irregular en¬ 
ergy bursts (see fig. [3] (f) for 200 <t< 350 and fig. [3] (j) for 
40 < t < 110) randomly occur on the relatively “smooth 
background” reminiscent of regimes shown on fig. 0(g) 
and (i). 

For 2073 < Ra < 2080, the family Bi gives rise to an 
attractor (B 2 ) displaying “in-out” intermittent behaviour 
[29] , i.e., the trajectory switches irregularly in time be¬ 
tween two main states, one is a time-periodic state corre¬ 
sponding to the destabilized A 2 , the other is chaotic, cor¬ 
responding to the destabilized Bi. Fig. [4] shows the time 
series, Poincare sections and velocity profiles for the at¬ 
tractor at Ra=2075; the trajectory in the phase space vis- 


p-3 





















R. Chertovskih et al. 



Fig. 3: Kinetic energy evolution and the Poincare section for 
the attractors of the family Bi for Ra=1880 (a), Ra=1910 (b), 
Ra=1990 (c), Ra=2000 (d), Ra=2010 (e), Ra=2035 (f), 
Ra=2040 (g), Ra=2050 (h), Ra=2060 (i), Ra=2070 (j). The 
axes are as in fig. [2] 


its the periodic state at 0 < t < 200 and 2100 <t< 3200. 



Fig. 4: Kinetic energy evolution, the Poincare sections and 
the vertival velocity profiles for the intermittent attractor at 
Ra=2075. On the top and (a) panels the same as on fig. [2] 
is shown. Points from (a) for the trajectories intersecting the 
Poincare plane for 400 < t < 1900 (b) and 2500 <t< 3000 (c). 
Four open circles in (c) show where the stable periodic orbit 
for Ra=2055 from branch A 2 intersects the Poincare plane. 
Isolines of V 3 on the horizontal mid-plane X 3 = 1/2 are shown 
for t = 1200 (d) and t = 2750 (e) (solid and dashed lines stand 
for positive and negative values, respectively). 


Thus, B 2 is formed by a crisis-like event involving the pre¬ 
viously destabilized A 2 and Bi. However, a comparison 
between the right panel of fig. 0 (j) and fig. I (a) reveals 
that B 2 is larger than the union of A 2 and Bi. The miss¬ 
ing component can be found from a set of initial condi¬ 
tions displaying chaotic transients for Ra < 2073. This 
set constitutes a chaotic saddle (CS), i.e., a nonattract¬ 


ing chaotic set responsible for chaotic transients 30 31 
Figs. §a) and (b) reveal the kinetic energy time series 
of two initial conditions at Ra=2070. Both trajectories 
exhibit chaotic behavior for about 100 time units, before 
they escape from the chaotic region toward the destabi¬ 
lized A 2 . Later, the same trajectories will converge to 
attractor Bi. The Poincare plot for the union of the tran¬ 
sient parts of both series is shown on fig. [ 5 ] (c) and resem¬ 
bles the Poincare plots of figs. [4] (b) and (c). From all this 
information, we conjecture the following scenario. Attrac¬ 
tor A 2 is destabilized in a boundary crisis at Ra~2055 af¬ 
ter a collision with CS, a chaotic saddle lying on the basin 
boundary between A 2 and Bi. The unstable set formed 
by the union of CS and the destabilized A 2 , then, collides 
with Bi at an interior crisis at Ra~2073 (IC1 in fig. [l]), 
leading to the intermittent attractor B 2 , where trajecto¬ 
ries alternate between phases where A 2 , CS and Bi can be 
identified. The highest energy bursts in fig. |4| correspond 
to CS. 

On increasing Ra, the intermittent attractor loses its 
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Fig. 5: Kinetic energy evolution of transients for Ra=2070 ((a) 
and (b)) and Poincare section for the transients (c). These 
transients are in the vicinity of a chaotic saddle in the phase 
space. The axes are as in fig. [2] 





Fig. 6 : Kinetic energy evolution for the attractors of the family 
B 3 for Ra=2140 (top panel) and Ra=2200 (bottom panel). 


stability, and trajectories for 2090 < Ra < 2130 are at¬ 
tracted to A 3 . However, chaotic transients reminiscent of 
the intermittent attractor B 2 remain all throughout this 
interval as a signature of a chaotic saddle in the back¬ 
ground of A 3 . For Ra > 2140, the stability of A 3 is lost, 
and a new chaotic attractor, B 3 , rises as (apparently) the 
sole attractor in the system. This attractor resembles B 2 
and the chaotic saddle (cf. the kinetic energy evolution for 
Ra=2075 from B 2 in fig. [4| and for Ra=2140 and Ra=2200 
from B 3 in fig. [ 6 ]). In the intermittent regimes of B 3 , the 
time spent near a state with regular behaviour is shorten¬ 
ing for increasing Ra (cf.regimes at Ra=2140, top panel 
of fig. [ 6 j and at Ra=2200, bottom panel of fig. [ 6 ]). Note 
that some of the regular phases of the intermittency are in 
the region previously occupied by A 3 . This corroborates 
another interior crisis scenario (IC2 in fig [l]) , whereby at¬ 
tractor A 3 collides with the background chaotic saddle to 
form an enlarged chaotic attractor B 3 , where intermittent 
switches between the former A 3 and the chaotic saddle 
take place. Such crisis-induced intermittency involving 
a destabilized attractor and a surrounding chaotic sad¬ 
dle was reported for a one-dimensional spatiotemporally 
chaotic system in [32] . 

The three largest Lyapunov exponents for the attractors 
from branches Bi and B 3 of fig. [T| are shown in fig. [7] as 
a function of Ra. Hyperchaos with at least three positive 
Lyapunov exponents is observed after the interior crisis of 
B 3 . Note the linear scaling of the exponents with Ra, a 
feature that has been observed in Rayleigh-Benard con¬ 
vection near the onset for the first exponent as a function 
of the reduced Rayleigh number by [33 . The gap between 


Fig. 7: The three largest Lyapunov exponents (vertical axis) 
for the attractors of the families Bi and B 3 . 


Bi and B 3 is filled by a chaotic saddle that evolves from 
the destabilized Bi. Although we have not computed the 
Lyapunov exponents of the chaotic saddle due to its re¬ 
pelling nature, we believe it should show a continuation of 
the linear scaling observed in fig. [7] 


Conclusions. — We have shown that the evolution 
from periodic convective states to hyperchaos in Rayleigh- 
Benard convection for P = 0.3 occurs through a sequence 
of local and global bifurcations as the Rayleigh number 
is increased, including three crises. First, a quasiperiodic 
attractor A 2 coexisting with another attractor Bi collides 
with a chaotic saddle (CS) in the basin boundary and loses 
stability due to a boundary crisis. Evidence for this is 
based on A 2 becoming transient and the disappearance 
of its basin of attraction. Then, the newly formed chaotic 
saddle (CS + destabilized A 2 ) collides with the chaotic at¬ 
tractor Bi in an interior crisis (IC1) leading to an enlarged 
attractor B 2 . The evidence is the intermittent switching 
between Bi, CS and A 2 . Next, attractor B 2 is destabilized 
and generates a chaotic saddle surrounding a quasiperiodic 
attractor A 3 (evidenced by B 2 becoming transient and the 
destruction of its basin of attraction). Finally, A 3 collides 
with the surrounding chaotic saddle in another interior cri¬ 
sis (IC2), which is evident from the intermittent switching 
between A 3 and CS in the time series. The enlarged at¬ 
tractor B 3 exhibits hyperchaos with at least three positive 
Lyapunov exponents. 

Although strictly speaking a crisis is a collision of a 
chaotic attractor with a saddle invariant set [34 , in a more 
general definition, a crisis can be described as a collision 
between any attractor and a saddle-type invariant set [35]. 
In the present letter, two of the identified crises involve the 
collision of a quasiperiodic attractor and a chaotic saddle. 

There is surprising similarity between our results and 


the route to chaos found in 


where two-dimensional 


Rayleigh-Benard convection was studied for P = 6 . 8 . 
However, the attractors found ibid, are of a limited in¬ 
terest in hydrodynamics, because two-dimensional convec¬ 
tive flows for the considered value of P are stable only in 
a narrow window of parameter values 


16 


. They are 

also not useful in the dynamo theory, since magnetic field 
generation by a planar flow is impossible by virtue of the 
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Zeldovich antidynamo theorem (see [l] for details). Note 
that our results also have some correlation with the bi¬ 
furcation diagram shown in 36 for zero Prandtl number 
near the onset of convection, where stationary, oscilating 
and chaotic regimes are found, although no hyperchaos 
was reported then, probably due to the low values of Ra. 
Intermittent and chaotic waves were also found in [37] for 
P = 0.2 as well as P = 0, with both free-slip and no-slip 
boundary conditions for small aspect ratio geometries. 

Experimental verification of some bifurcations reported 
here have been found in the past for Rayleigh-Benard con¬ 
vection, such as abrupt transitions from quasiperiodic to 
chaotic attractors [38], a route to chaos via intermittency 
[39] , and quasiperiodic attractors with three basic frequen¬ 
cies [40] , 

Regarding future works, note that we have not fully 
characterised all bifurcations presented and, in particular, 
the destabilisation of B 2 (supposedly another boundary 
crisis) was not investigated. Additionally, we do not claim 
to have found all the attractors present in the range of 
Ra studied. Continuing the line of research started in 
we plan to add rotation and magnetic field, 
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whereby one can study mechanisms of the convective elec¬ 
tromagnetic processes in the liquid outer core of the Earth. 
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